<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><script type="text/javascript" src="https://cdn.jsdelivr.net/gh/opencobra/cobratoolbox@ffa0229fc0c01c9236bb7e961f65712443277719/latest/_static/js/iframeResizer.contentWindow.min.js"></script><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Adding biological constraints to a flux balance model</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: normal; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 15px 10px 5px 4px; padding: 0px; line-height: 18px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 17px; font-weight: bold; text-align: left;  }
.S2 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S3 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S4 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: normal; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S5 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: center;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S7 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: normal;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: normal;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: normal;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}
.S8 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S13 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S14 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S15 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: normal;  }
.variableValue { width: 100% !important; }
.embeddedOutputsMatrixElement,.eoOutputWrapper .matrixElement {    min-height: 18px;    box-sizing: border-box;}
.embeddedOutputsMatrixElement .matrixElement,.eoOutputWrapper  .matrixElement,.rtcDataTipElement .matrixElement {    position: relative;}
.matrixElement .variableValue,.rtcDataTipElement .matrixElement .variableValue {    white-space: normal;    display: inline-block;    vertical-align: top;    overflow: hidden;}
.embeddedOutputsMatrixElement.inlineElement {}
.embeddedOutputsMatrixElement.inlineElement .topHeaderWrapper {    display: none;}
.embeddedOutputsMatrixElement.inlineElement .veTable .body {    padding-top: 0 !important;    max-height: 100px;}
.inlineElement .matrixElement {    max-height: 300px;}
.embeddedOutputsMatrixElement.rightPaneElement {}
.rightPaneElement .matrixElement,.rtcDataTipElement .matrixElement {    overflow: hidden;    padding-left: 9px;}
.rightPaneElement .matrixElement {    margin-bottom: -1px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer,.eoOutputWrapper .matrixElement .valueContainer,.rtcDataTipElement .matrixElement .valueContainer {    white-space: nowrap;    margin-bottom: 3px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis.hide,.embeddedOutputsMatrixElement .matrixElement .verticalEllipsis.hide,.eoOutputWrapper .matrixElement .valueContainer .horizontalEllipsis.hide,.eoOutputWrapper .matrixElement .verticalEllipsis.hide,.rtcDataTipElement .matrixElement .valueContainer .horizontalEllipsis.hide,.rtcDataTipElement .matrixElement .verticalEllipsis.hide {    display: none;}
.embeddedOutputsVariableMatrixElement .matrixElement .valueContainer.hideEllipses .verticalEllipsis, .embeddedOutputsVariableMatrixElement .matrixElement .valueContainer.hideEllipses .horizontalEllipsis {    display:none;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis,.eoOutputWrapper .matrixElement .valueContainer .horizontalEllipsis {    margin-bottom: -3px;}
.eoOutputWrapper .embeddedOutputsVariableMatrixElement .matrixElement .valueContainer {    cursor: default !important;}
/* * Ellipses as base64 for HTML export. */.matrixElement .horizontalEllipsis,.rtcDataTipElement .matrixElement .horizontalEllipsis {    display: inline-block;    margin-top: 3px;    /* base64 encoded version of images-liveeditor/HEllipsis.png */    width: 30px;    height: 12px;    background-repeat: no-repeat;    background-image: url("");}
.matrixElement .verticalEllipsis,.textElement .verticalEllipsis,.rtcDataTipElement .matrixElement .verticalEllipsis,.rtcDataTipElement .textElement .verticalEllipsis {    margin-left: 35px;    /* base64 encoded version of images-liveeditor/VEllipsis.png */    width: 12px;    height: 30px;    background-repeat: no-repeat;    background-image: url("");}
.S16 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }</style></head><body><div class = rtcContent><h1  class = 'S0'><span>Adding biological constraints to a flux balance model</span></h1><h3  class = 'S1'><span style=' font-weight: bold;'>Note: This tutorial is a draft and needs completion. Contributions welcome!</span></h3><h2  class = 'S2'><span>Authors: Diana C. El Assal and Ronan M.T. Fleming, </span><span style=' font-weight: bold;'> </span><span>Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Luxembourg.</span></h2><h2  class = 'S2'><span>Reviewers: </span></h2><div  class = 'S3'><span>Anne Richelle, Systems Biology and Cell engineering, University of California San Diego</span></div><div  class = 'S3'><span>Almut Heinken, Molecular Systems Physiology Group, LCSB, University of Luxembourg</span></div><h2  class = 'S4'><span>INTRODUCTION</span></h2><div  class = 'S3'><span>A metabolic model can be converted into a condition-specific model based on the imposition of experimentally derived constraints. Constraints can be defined, for example, by imposing upper and lower flux bounds for each reaction. There are several types of constraints that can be imposed in a metabolic model and that represent specific intra- and extracellular conditions, such as biomass maintenance requirements, environmental constraints, or maximum enzyme capacities. </span></div><div  class = 'S3'><span>In general, biomass constraints [1] are added as part of a biomass reaction. In some instances, however, a cell-type (e.g. neurons) does not divide, but is only required to turn over its biomass components. This tutorial is particularly relevant for such cases. Turnover rates are commonly expressed as half-lives (</span><span texencoding="{t_{1/2}" style="vertical-align:-6px"><img src="" width="22" height="21" /></span><span>) and represent the time required for half of the biomass precursor to be replaced [2]. </span></div><div  class = 'S3'><span>Using the experimental literature, metabolite </span><span texencoding="{t_{1/2}" style="vertical-align:-6px"><img src="" width="22" height="21" /></span><span> were collected and converted into turnover rates (</span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">λ</span><span>):</span></div><div  class = 'S5'><span texencoding="\text{(1)}\ \lambda=\frac{ln(2)}{t_{1/2}}\\
" style="vertical-align:-24px"><img src="" width="84.5" height="59" /></span></div><div  class = 'S3'><span>We consider a biochemical network of </span><span> </span><span>m</span><span> </span><span> molecular species and </span><span> </span><span>n</span><span> </span><span> biochemical reactions. The biochemical network is mathematically represented by a stoichiometric matrix </span><span texencoding="S\in\mathcal{Z}^{m\times n}" style="vertical-align:-5px"><img src="" width="63.5" height="19" /></span><span>.  After calculating </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">λ</span><span>, we integrate it into the steady-state equation:</span></div><div  class = 'S5'><span texencoding="\text{(2)}\ Sv=\frac{dx}{dt}\\" style="vertical-align:-22px"><img src="" width="74.5" height="56" /></span></div><div  class = 'S3'><span>The steady-state flux vector </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">v</span><span> and the change in abundance over time (</span><span texencoding="\frac{dx}{dt}\\
" style="vertical-align:-22px"><img src="" width="17" height="56" /></span><span>) share the same units, with  </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">x</span><span>  being the abundance of the corresponding biomass precursor. Equation (2) can thus be re-written: </span></div><div  class = 'S5'><span texencoding="
\text{(3)}\ Sv=\frac{dx}{dt}=\lambda x\\
" style="vertical-align:-22px"><img src="" width="107.5" height="56" /></span></div><h2  class = 'S4'><span>PROCEDURE</span></h2><div  class = 'S3'><span>Initialize the Cobra Toolbox using the </span><span style=' font-family: monospace;'>initCobraToolbox</span><span> function.</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span >initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% false, as we don't want to update</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="4BE0BE33" data-testid="output_0" data-width="428" data-height="885" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">      _____   _____   _____   _____     _____     |
     /  ___| /  _  \ |  _  \ |  _  \   / ___ \    |   COnstraint-Based Reconstruction and Analysis
     | |     | | | | | |_| | | |_| |  | |___| |   |   The COBRA Toolbox - 2017
     | |     | | | | |  _  { |  _  /  |  ___  |   |
     | |___  | |_| | | |_| | | | \ \  | |   | |   |   Documentation:
     \_____| \_____/ |_____/ |_|  \_\ |_|   |_|   |   <a href="http://opencobra.github.io/cobratoolbox" style="white-space: normal; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">http://opencobra.github.io/cobratoolbox</a>
                                                  | 

 &gt; Checking if git is installed ...  Done.
 &gt; Checking if the repository is tracked using git ...  Done.
 &gt; Checking if curl is installed ...  Done.
 &gt; Checking if remote can be reached ...  Done.
 &gt; Initializing and updating submodules ... Done.
 &gt; Adding all the files of The COBRA Toolbox ...  Done.
 &gt; Define CB map output... set to svg.
 &gt; Retrieving models ...   Done.
 &gt; TranslateSBML is installed and working properly.
 &gt; Configuring solver environment variables ...
   - [*---] ILOG_CPLEX_PATH: C:\Program Files\IBM\ILOG\CPLEX_Studio1271\cplex\matlab\x64_win64
   - [*---] GUROBI_PATH: C:\gurobi650\win64\matlab
   - [*---] TOMLAB_PATH: C:\tomlab\
   - [----] MOSEK_PATH :  --&gt; set this path manually after installing the solver ( see <a href="https://opencobra.github.io/cobratoolbox/docs/solvers.html" style="white-space: normal; font-style: normal; color: rgb(0, 95, 206); font-size: 12px;">instructions</a> )
   Done.
 &gt; Checking available solvers and solver interfaces ... Done.
 &gt; Setting default solvers ... Done.
 &gt; Saving the MATLAB path ... Done.
   - The MATLAB path was saved in the default location.

 &gt; Summary of available solvers and solver interfaces

					Support           LP 	 MILP 	   QP 	 MIQP 	  NLP
	----------------------------------------------------------------------
	cplex_direct 	active        	    0 	    0 	    0 	    0 	    -
	dqqMinos     	active        	    0 	    - 	    - 	    - 	    -
	glpk         	active        	    1 	    1 	    - 	    - 	    -
	gurobi       	active        	    1 	    1 	    1 	    1 	    -
	ibm_cplex    	active        	    1 	    1 	    1 	    - 	    -
	matlab       	active        	    1 	    - 	    - 	    - 	    1
	mosek        	active        	    0 	    0 	    0 	    - 	    -
	pdco         	active        	    1 	    - 	    1 	    - 	    -
	quadMinos    	active        	    0 	    - 	    - 	    - 	    0
	tomlab_cplex 	active        	    1 	    1 	    1 	    1 	    -
	qpng         	passive       	    - 	    - 	    1 	    - 	    -
	tomlab_snopt 	passive       	    - 	    - 	    - 	    - 	    1
	gurobi_mex   	legacy        	    0 	    0 	    0 	    0 	    -
	lindo_old    	legacy        	    0 	    - 	    - 	    - 	    -
	lindo_legacy 	legacy        	    0 	    - 	    - 	    - 	    -
	lp_solve     	legacy        	    1 	    - 	    - 	    - 	    -
	opti         	legacy        	    0 	    0 	    0 	    0 	    0
	----------------------------------------------------------------------
	Total        	-             	    7 	    4 	    5 	    2 	    2

 + Legend: - = not applicable, 0 = solver not compatible or not installed, 1 = solver installed.


 &gt; You can solve LP problems using: 'glpk' - 'gurobi' - 'ibm_cplex' - 'matlab' - 'pdco' - 'tomlab_cplex' - 'lp_solve' 
 &gt; You can solve MILP problems using: 'glpk' - 'gurobi' - 'ibm_cplex' - 'tomlab_cplex' 
 &gt; You can solve QP problems using: 'gurobi' - 'ibm_cplex' - 'pdco' - 'tomlab_cplex' - 'qpng' 
 &gt; You can solve MIQP problems using: 'gurobi' - 'tomlab_cplex' 
 &gt; You can solve NLP problems using: 'matlab' - 'tomlab_snopt' 

 &gt; Checking for available updates ...
 &gt; The COBRA Toolbox is up-to-date.</div></div></div></div></div><h2  class = 'S4'><span style=' font-weight: bold;'>Setting the </span><span>optimization</span><span style=' font-weight: bold;'> solver</span></h2><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span >changeCobraSolver(</span><span style="color: rgb(170, 4, 249);">'gurobi'</span><span >,</span><span style="color: rgb(170, 4, 249);">'LP'</span><span >);</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="9F83CB76" data-testid="output_1" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> &gt; Gurobi interface added to MATLAB path.</div></div></div></div></div><div  class = 'S8'><span>Here, we use Recon2.0 model (distributed by the toolbox) for illustration, although any model can be used. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >modelFileName = </span><span style="color: rgb(170, 4, 249);">'Recon2.0model.mat'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >modelDirectory = getDistributedModelFolder(modelFileName); </span><span style="color: rgb(2, 128, 9);">%Look up the folder for the distributed Models.</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >modelFileName= [modelDirectory filesep modelFileName]; </span><span style="color: rgb(2, 128, 9);">% Get the full path. Necessary to be sure, that the right model is loaded</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >model = readCbModel(modelFileName);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >modelOrig = model;</span></span></div></div></div><h2  class = 'S4'><span style=' font-weight: bold;'>1. Environmental constraints</span></h2><div  class = 'S5'><img class = "imageNode" src = "" width = "400" height = "98" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S3'><span>Environmental constraints are typically related to nutrient availability (e.g., glucose and oxygen). They can be defined using the function </span><span style=' font-style: italic;'>changeRxnBounds</span><span> to set the minimal and maximal uptake and/or secretion rates possible in a specific condition. For example, in the caudate-putamen of the conscious rat, glucose consumption rate was found to range between -12.00</span><span> and </span><span>-11.58</span><span>  </span><span texencoding="\mu mol/gDW/hr" style="vertical-align:-5px"><img src="" width="95.5" height="19" /></span><span>[3]. Therefore, the lower bound of the glucose exchange reaction (</span><a href = "http://www.vmh.life/#reaction/EX_glc(e)"><span>EX_glc(e)</span></a><span>) can be set as follows:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >modelConstrained = model;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >modelConstrained.c = 0*modelConstrained.c; </span><span style="color: rgb(2, 128, 9);">% remove any objective function</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >modelConstrained = changeRxnBounds(modelConstrained, </span><span style="color: rgb(170, 4, 249);">'EX_glc(e)'</span><span >, -12, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><div  class = 'S8'><span>Optionally, to further constrain the model, an upper bound can also be imposed to force the model to take up between 11.58 and 12 units of glucose </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >modelConstrained = changeRxnBounds(modelConstrained, </span><span style="color: rgb(170, 4, 249);">'EX_glc(e)'</span><span >, -11.58, </span><span style="color: rgb(170, 4, 249);">'u'</span><span >);</span></span></div></div></div><h2  class = 'S4'><span style=' font-weight: bold;'>2. Internal enzymatic constraints </span></h2><div  class = 'S5'><img class = "imageNode" src = "" width = "400" height = "99" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S3'><span>By convention, the bounds set on reaction rates in a metabolic model range from -1000 to 1000 and from 0 to 1000 for reversible and irreversible reactions, respectively [4]. Actually, the rate of a reaction is related to the activity of the enzyme catalyzing this reaction.  Therefore, internal enzymatic constraints can be used to define the maximum capacity of a specified enzyme to catalyze a reaction (</span><span style=' font-style: italic;'>Vmax</span><span>). For example, assuming that the reaction catalyzed by fructose-bisphosphate aldolase</span><span> (</span><a href = "http://www.vmh.life/#reaction/FBA"><span>FBA</span></a><span>) has a </span><span style=' font-style: italic;'>Vmax </span><span>of 128 units in our specific cell type, we can then add this constraint on the corresponding internal reaction FBA, as an upper bound. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >modelConstrained = changeRxnBounds(modelConstrained, </span><span style="color: rgb(170, 4, 249);">'FBA'</span><span >, 128, </span><span style="color: rgb(170, 4, 249);">'u'</span><span >);</span></span></div></div></div><div  class = 'S8'><span>Optionally, if the reaction is reversible, the same constraint can be set as the lower bound, but with opposite signs.</span></div><div  class = 'S5'><img class = "imageNode" src = "" width = "400" height = "97" alt = "" style = "vertical-align: baseline"></img></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >modelConstrained = changeRxnBounds(modelConstrained, </span><span style="color: rgb(170, 4, 249);">'FBA'</span><span >, -128, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><div  class = 'S8'><span></span></div><h2  class = 'S4'><span style=' font-weight: bold;'>3. Constraints associated with biomass</span></h2><div  class = 'S3'><span>In general, biomass constraints [1] are added as part of a biomass reaction by defining stoichiometric coefficients for each biomass precursor. For dividing cell types, the generic human biomass reaction available in Recon2 is formulated as follows: </span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span >printRxnFormula(modelConstrained, </span><span style="color: rgb(170, 4, 249);">'biomass_reaction'</span><span >);</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="C70B053F" data-testid="output_2" data-width="428" data-height="18" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">biomass_reaction	20.6508 h2o[c] + 20.7045 atp[c] + 0.385872 glu_L[c] + 0.352607 asp_L[c] + 0.036117 gtp[c] + 0.279425 asn_L[c] + 0.505626 ala_L[c] + 0.046571 cys_L[c] + 0.325996 gln_L[c] + 0.538891 gly[c] + 0.392525 ser_L[c] + 0.31269 thr_L[c] + 0.592114 lys_L[c] + 0.35926 arg_L[c] + 0.153018 met_L[c] + 0.023315 pail_hs[c] + 0.039036 ctp[c] + 0.154463 pchol_hs[c] + 0.055374 pe_hs[c] + 0.020401 chsterol[c] + 0.002914 pglyc_hs[c] + 0.011658 clpn_hs[c] + 0.053446 utp[c] + 0.009898 dgtp[n] + 0.009442 dctp[n] + 0.013183 datp[n] + 0.013091 dttp[n] + 0.275194 g6p[c] + 0.126406 his_L[c] + 0.159671 tyr_L[c] + 0.286078 ile_L[c] + 0.545544 leu_L[c] + 0.013306 trp_L[c] + 0.259466 phe_L[c] + 0.412484 pro_L[c] + 0.005829 ps_hs[c] + 0.017486 sphmyln_hs[c] + 0.352607 val_L[c] 	-&gt;	20.6508 adp[c] + 20.6508 h[c] + 20.6508 pi[c] </div></div></div></div></div><h2  class = 'S2'><span style=' font-weight: bold;'>3.1 Biomass reaction</span></h2><div  class = 'S3'><span>20.6508 h2o[c] + 20.7045 atp[c] + 0.38587 glu_L[c] + 0.35261 asp_L[c] + 0.036117 gtp[c] + 0.50563 ala_L[c] + 0.27942 asn_L[c] + 0.046571 cys_L[c] + 0.326 gln_L[c] + 0.53889 gly[c] + 0.39253 ser_L[c] + 0.31269 thr_L[c] + 0.59211 lys_L[c] + 0.35926 arg_L[c] + 0.15302 met_L[c] + 0.023315 pail_hs[c] + 0.039036 ctp[c] + 0.15446 pchol_hs[c] + 0.055374 pe_hs[c] + 0.020401 chsterol[c] + 0.002914 pglyc_hs[c] + 0.011658 clpn_hs[c] + 0.009898 dgtp[n] + 0.009442 dctp[n] + 0.013183 datp[n] + 0.053446 utp[c] + 0.013091 dttp[n] + 0.27519 g6p[c] + 0.12641 his_L[c] + 0.15967 tyr_L[c] + 0.28608 ile_L[c] + 0.54554 leu_L[c] + 0.013306 trp_L[c] + 0.25947 phe_L[c] + 0.41248 pro_L[c] + 0.005829 ps_hs[c] + 0.017486 sphmyln_hs[c] + 0.35261 val_L[c]  -&gt; 20.6508 h[c] + 20.6508 adp[c] + 20.6508 pi[c] </span></div><div  class = 'S3'><span>Any changes or adaptations can be introduced by adding a new formulation of the biomass function, using the function </span><span style=' font-style: italic;'>addReaction</span><span>. For example, one can add the following new biomass reaction named </span><span style=' font-style: italic;'>biomassReactionLipids</span><span>:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span >modelConstrained = addReaction(modelConstrained, </span><span style="color: rgb(170, 4, 249);">'biomasReactionLipids'</span><span >,  </span><span style="color: rgb(170, 4, 249);">'20.6508 h2o[c] + 20.7045 atp[c] + 0.15446 pchol_hs[c] + 0.055374 pe_hs[c] + 0.020401 chsterol[c] + 0.011658 clpn_hs[c] + 0.005829 ps_hs[c] + 0.017486 sphmyln_hs[c] -&gt; 20.6508 h[c] + 20.6508 adp[c] + 20.6508 pi[c]'</span><span >);</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="753D976E" data-testid="output_3" data-width="428" data-height="31" data-hashorizontaloverflow="true" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">test
biomasReactionLipids	20.6508 h2o[c] + 20.7045 atp[c] + 0.15446 pchol_hs[c] + 0.055374 pe_hs[c] + 0.020401 chsterol[c] + 0.011658 clpn_hs[c] + 0.005829 ps_hs[c] + 0.017486 sphmyln_hs[c] 	-&gt;	20.6508 adp[c] + 20.6508 h[c] + 20.6508 pi[c] </div></div></div></div></div><h2  class = 'S4'><span style=' font-weight: bold;'>4. Biomass maintenance constraints</span></h2><div  class = 'S3'><span>To represent biomass maintenance (e.g. in neurons), the minimal biomass maintenance requirements can be used to set the corresponding constraints. Using the neurobiochemical literature, the degradation pathways for each biomass precursor has to be identified and the corresponding first reactions of these degradation pathways need to be mapped to ReconX. Using the fractional composition and the turnover rate of each biomass precursor, corresponding reaction rates (</span><span texencoding="\mu mol/gDW/hr" style="vertical-align:-5px"><img src="" width="95.5" height="19" /></span><span>) are calculated as described above and in Table 1. These reaction rates represent the minimal requirements for biomass maintenance of neurons in the human grey matter and must therefore be imposed as a lower bound on the corresponding degradation reaction(s) of the different lipids, amino acids, and nucleic acids. </span></div><div  class = 'S3'><img class = "imageNode" src = "" width = "400" height = "94" alt = "" style = "vertical-align: baseline"></img><a href = "#null"><span></span></a></div><div  class = 'S3'><span style=' font-weight: bold;'>Table 1: The minimum metabolic maintenance requirement for neurons.</span><span> </span><span>This is a coarse-grained approximation of neuronal lipid, amino acid, and nucleic acid maintenance requirements converted into </span><span texencoding="\mu mol/gDW/hr" style="vertical-align:-5px"><img src="" width="95.5" height="19" /></span><span>.  </span></div><div  class = 'S5'><img class = "imageNode" src = "" width = "614" height = "417" alt = "" style = "vertical-align: baseline"></img></div><div  class = 'S3'><span>**cardiolipin is also known as diphosphatidylglycerol</span></div><div  class = 'S3'><span style=' font-weight: bold; font-style: italic;'>Calculation example of the minimal cholesterol maintenance requirement </span></div><div  class = 'S3'><span style=' font-weight: bold; font-style: italic;'>i. Identify metabolite abundance</span></div><div  class = 'S3'><span>Assuming that the specific tissue type has a total dry weight lipid composition of 39.6%. This means that there is 0.396gLipid/gDWtissue. If cholesterol has a molar composition of 31.3%, then in total there is 0.124g cholesterol per gDW of tissue.</span></div><div  class = 'S5'><span texencoding="\frac{31.3*39.6*1gDW}{100*100}=0.124gDW\\\
" style="vertical-align:-22px"><img src="" width="209.5" height="56" /></span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >abundance = (31.3*39.6*1)/(100*100);</span></span></div></div></div><div  class = 'S8'><span style=' font-weight: bold; font-style: italic;'>ii. Calculate the molar abundance</span></div><div  class = 'S3'><span>In the experimental literature, cholesterol was also found to have a molar mass (</span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">M</span><span>) of 386</span><span texencoding="g/mol" style="vertical-align:-5px"><img src="" width="40.5" height="19" /></span><span>. Using equation (4), we can now convert the abundance (</span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">m</span><span>) into molar units (</span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">n</span><span>). </span></div><div  class = 'S5'><span texencoding="
\text{(4)}\ n=\frac{m}{M}
" style="vertical-align:-15px"><img src="" width="68.5" height="35" /></span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >M = 386; </span><span style="color: rgb(2, 128, 9);">%g/mol</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >n = (abundance*1000000)/M; </span><span style="color: rgb(2, 128, 9);">%micromol</span></span></div></div></div><div  class = 'S8'><span style=' font-weight: bold; font-style: italic;'>iii. Calculate the corresponding flux value</span></div><div  class = 'S3'><span>Finally, we know that in the brain, cholesterol has a very slow turnover and a </span><span texencoding="{t_{1/2}" style="vertical-align:-6px"><img src="" width="22" height="21" /></span><span> of 4320 hours. Using equation (3), we can now calculate the minimal cholesterol maintenance requirement in flux units (</span><span texencoding="v_1" style="vertical-align:-6px"><img src="" width="14.5" height="20" /></span><span>). </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >halfLife = 4320;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >turnover = log(2)/halfLife;</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S13'><span style="white-space: normal"><span >v1 = n * turnover</span></span></div><div  class = 'S7'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>v1 = 0.0515</div></div></div></div><div  class = 'S8'><span>The minimal cholesterol maintenance requirement was calculated to be 0.0515</span><span texencoding="\mu mol/gDW/hr" style="vertical-align:-5px"><img src="" width="95.5" height="19" /></span><span>(Table 1). This value can now be used as a lower bound in the corresponding reaction.</span></div><h2  class = 'S2'><span style=' font-weight: bold;'>4.1. </span><span>Identification of degradation reactions for a biomass maintenance precursor</span></h2><div  class = 'S3'><span>As previously mentioned, the degradation pathways for each biomass precursor can be identified using literature and the minimal maintenance requirements defined in section 4. can be used to constrain the first reactions of these degradation pathways. However, the identification of these reactions and the set-up of the associated constraints is not always straightforward. The following section presents the different common cases that can be encountered.</span></div><h2  class = 'S2'><span>A. Single irreversible degradation reaction</span></h2><div  class = 'S3'><span>In cases where only a single irreversible degradation reaction exists for a biomass maintenance precursor, the imposition of the constraint is straightforward. For example, the major cholesterol excretion pathway in the brain involves the hydroxylation of cholesterol into the oxysterol 24-hydroxycholesterol. Only a subset of neurons express this 24-hydroxylase enzyme (</span><a href = "http://www.vmh.life/#human/all/P45046A1r"><span>P45046A1r</span></a><span>) and it is mainly found in dendrites and somata, rather than in axons or presynaptic terminals (reviewed in [5]). </span></div><div  class = 'S3'><span>Therefore, add a lower bound (</span><span texencoding="v_1" style="vertical-align:-6px"><img src="" width="14.5" height="20" /></span><span>, calculated above) on P45046A1r (Table 1). </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >modelConstrained = changeRxnBounds(modelConstrained, </span><span style="color: rgb(170, 4, 249);">'P45046A1r'</span><span >, v1, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><h2  class = 'S2'><span>B. Single degradation reaction does not exist biochemically</span></h2><div  class = 'S3'><span>A degradation reaction might not exist for a given biomass maintenance precursor. For example, the phospholipid cardiolipin is mainly present in the inner mitochondrial membrane, where it regulates the stability of the mitochondrial membrane protein complexes [6]. As part of mitochondria, cardiolipin reaches the lysosome during macroautophagy (reviewed in [7]). It is then degraded to form the negatively charged bis(monoacylglycero)phosphate (BMP) on internal membranes.</span><span>  </span></div><div  class = 'S3'><span>If the corresponding demand reaction does not exist in the model, a demand reaction can be added using the function </span><span style=' font-style: italic;'>addDemandReaction</span><span>. </span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span >modelConstrained = addDemandReaction(modelConstrained, </span><span style="color: rgb(170, 4, 249);">'clpn_hs[c]'</span><span >); </span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement" uid="12CF9B18" data-testid="output_5" data-width="428" data-height="18" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">DM_clpn_hs[c]	clpn_hs[c] 	-&gt;	</div></div></div></div></div><div  class = 'S8'><span>Now, the constraint for cardiolipin (Table 1) can be imposed on the corresponding demand reaction</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >modelConstrained = changeRxnBounds(modelConstrained, </span><span style="color: rgb(170, 4, 249);">'DM_clpn_hs[c]'</span><span >, 0.001, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div></div><h2  class = 'S2'><span>C. Single reversible degradation reaction</span></h2><div  class = 'S3'><span>In cases where the biomass precursor degradation reaction is reversible, it first needs to be split into two irreversible reactions. To this end, define the set of reversible degradation reactions (sRxns) and convert them into two irreversible reactions (i.e., </span><span style=' font-style: italic;'>sRXN_b</span><span> and </span><span style=' font-style: italic;'>sRXN_f</span><span>, respectively backward and forward reactions) using </span><span style=' font-style: italic;'>convertToIrreversible:</span></div><div  class = 'S3'><span style=' font-weight: bold;'>i. Split the reversible reactions into irreversible</span></div><div  class = 'S3'><span>Select a set of reversible degradation reactions </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >sRxns = {</span><span style="color: rgb(170, 4, 249);">'ASPTA'</span><span >; </span><span style="color: rgb(170, 4, 249);">'GHMT2r'</span><span >};</span></span></div></div></div><div  class = 'S8'><span>Copy the original model, except split a specific list of reversible degradation reactions into irreversible. </span></div><div  class = 'S3'><span>Split sRxns into irreversible reactions</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >[modelIrrev] = convertToIrreversible(modelConstrained,</span><span style="color: rgb(170, 4, 249);">'sRxns'</span><span >,sRxns);</span></span></div></div></div><div  class = 'S8'><span>You can check if the conversion as been done properly by searching the split reactions</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >j=1:length(sRxns)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >isempty(findRxnIDs(modelIrrev, [sRxns{j} </span><span style="color: rgb(170, 4, 249);">'_f'</span><span >]))</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >        error(</span><span style="color: rgb(170, 4, 249);">'Forward reaction not found'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >isempty(findRxnIDs(modelIrrev, [sRxns{j} </span><span style="color: rgb(170, 4, 249);">'_b'</span><span >]))</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >        error(</span><span style="color: rgb(170, 4, 249);">'Reverse reaction not found'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S8'><span style=' font-weight: bold;'>ii. Impose the calculated constraints</span></div><div  class = 'S3'><span>Examples are given for aspartate and glutamate (Table 1). </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >constraints = [1.590; 1.146];</span></span></div></div></div><div  class = 'S8'><span>You can also identify the new reaction names as follows:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span >rxns = setdiff(modelIrrev.rxns, modelConstrained.rxns)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="D9ED7F56" data-testid="output_6" data-width="428" data-height="76" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">rxns = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'ASPTA_b'
    'ASPTA_f'
    'GHMT2r_b'
    'GHMT2r_f'
</div></div></div></div></div></div><div  class = 'S8'><span>Using this list (rxns), manually identify the corresponding reactions that should be constrained</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >splitRxns= {</span><span style="color: rgb(170, 4, 249);">'ASPTA_f'</span><span >; </span><span style="color: rgb(170, 4, 249);">'GHMT2r_f'</span><span >};</span></span></div></div></div><div  class = 'S8'><span> Identify the indices of these split reactions in the new model, using </span><span style=' font-style: italic;'>findRxnIDs</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >ind = findRxnIDs(modelIrrev, splitRxns);</span></span></div></div></div><div  class = 'S8'><span>Now, impose the constraints using a for loop. Note that you can easily account for experimental errors by defining a percentage error (e.g., </span><span style=' font-style: italic;'>expError = 0.25</span><span>) for the constraint values. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >expError = 0.25;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >modelConstrained = modelIrrev;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >i = 1:length(splitRxns)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >    modelConstrained = changeRxnBounds(modelConstrained, splitRxns{i,1},</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >        constraints(i,1)-constraints(i,1)*expError, </span><span style="color: rgb(170, 4, 249);">'l'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><h2  class = 'S2'><span>D. Multiple irreversible degradation reactions</span></h2><div  class = 'S3'><span>In some cases, several degradation pathways may be available for one biomass precursor. For example, in the brain, phosphatidylcholine (PC) can be degraded by 3 different metabolic pathways [8]:</span></div><ul  class = 'S14'><li  class = 'S15'><a href = "http://www.vmh.life/#reaction/PCHOLP_hs"><span>PCHOLP_hs</span></a><span>: Phospholipase D acts on the choline/phosphate bond of PC to form choline and phosphatidic acid.</span></li><li  class = 'S15'><a href = "http://www.vmh.life/#reaction/PLA2_2"><span>PLA2_2</span></a><span>: Phospholipase A2 acts on the bond between the fatty acid and the hydroxyl group of PC to form a fatty acid (e.g. arachidonic acid or docosahexaenoic acid) and lysophosphatidylcholine.</span><span> </span></li><li  class = 'S15'><a href = "http://www.vmh.life/#reaction/SMS"><span>SMS</span></a><span>: Ceramide and PC can also be converted to sphingomyelin by sphingomyelin synthetase.</span><span> </span></li></ul><div  class = 'S3'><span>Define the set of potential reactions associated with the degradation of PC</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >multipleRxnList={</span><span style="color: rgb(170, 4, 249);">'PCHOLP_hs'</span><span >, </span><span style="color: rgb(170, 4, 249);">'PLA2_2'</span><span >, </span><span style="color: rgb(170, 4, 249);">'SMS'</span><span >};</span></span></div></div></div><div  class = 'S8'><span>Make sure that all the reactions are irreversible. The lower bounds should be 0 and the upper bounds 1000.  </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >modelConstrained.lb(findRxnIDs(modelConstrained, multipleRxnList));</span></span></div></div></div><div  class = 'S8'><span>ans =</span></div><div  class = 'S3'><span>     0</span></div><div  class = 'S3'><span>     0</span></div><div  class = 'S3'><span>     0</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S12'><span style="white-space: normal"><span >modelConstrained.ub(findRxnIDs(modelConstrained, multipleRxnList));</span></span></div></div></div><div  class = 'S8'><span>ans =</span></div><div  class = 'S3'><span>        1000</span></div><div  class = 'S3'><span>        1000</span></div><div  class = 'S3'><span>        1000</span></div><div  class = 'S3'><span>Constrain the weighted sum of fluxes to be above a lower bound (e.g. value of the maintenance requirement of PC in Table 1: d = 2.674 </span><span texencoding="\mu mol/gDW/hr" style="vertical-align:-5px"><img src="" width="95.5" height="19" /></span><span>). The weight for each reaction are defined in </span><span style=' font-style: italic;'>c</span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >c=[1,1,1];</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >d=2.674; </span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >ineqSense=</span><span style="color: rgb(170, 4, 249);">'G'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >modelConstrainedAb=constrainRxnListAboveBound(modelConstrained,multipleRxnList,c,d,ineqSense);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >rxnInd = findRxnIDs(modelConstrainedAb, multipleRxnList);</span></span></div></div></div><div  class = 'S8'><span>Check the constraints are there</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >[nMet,nRxn]=size(modelConstrainedAb.S);</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S13'><span style="white-space: normal"><span >modelConstrainedAb</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement scrollableOutput" uid="B08D2DC2" data-testid="output_7" data-width="428" data-height="608" data-hashorizontaloverflow="false" style="width: 458px; max-height: 261px; white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">modelConstrainedAb = </span></div><div style="white-space: normal; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">                      S: [5063×7444 double]
                   rxns: {7444×1 cell}
                     lb: [7444×1 double]
                     ub: [7444×1 double]
                    rev: [7440×1 double]
                      c: [7444×1 double]
             rxnGeneMat: [7444×2194 double]
                  rules: {7444×1 cell}
                  genes: {2194×1 cell}
                grRules: {7444×1 cell}
             subSystems: {7444×1 cell}
               rxnNames: {7444×1 cell}
              rxnKeggID: {7440×1 cell}
    rxnConfidenceEcoIDA: {7440×1 cell}
    rxnConfidenceScores: {7440×1 cell}
             rxnsboTerm: {7440×1 cell}
          rxnReferences: {7440×1 cell}
           rxnECNumbers: {7440×1 cell}
               rxnNotes: {7440×1 cell}
                   mets: {5063×1 cell}
                      b: [5063×1 double]
               metNames: {5063×1 cell}
            metFormulas: {5063×1 cell}
              metCharge: [5063×1 double]
             metCHEBIID: {5063×1 cell}
              metKeggID: {5063×1 cell}
           metPubChemID: {5063×1 cell}
         metInchiString: {5063×1 cell}
         metHepatoNetID: {5063×1 cell}
              metEHMNID: {5063×1 cell}
            ExchRxnBool: [7440×1 logical]
              EXRxnBool: [7440×1 logical]
              DMRxnBool: [7440×1 logical]
            SinkRxnBool: [7440×1 logical]
            SIntRxnBool: [7440×1 logical]
                metHMDB: {5063×1 cell}
                modelID: 'Recon2.0model'
                 csense: [5064×1 char]
                  match: [7444×1 double]
        reversibleModel: 0
                      C: [1×7444 double]
                      d: 2.6740
</div></div></div></div></div></div><div  class = 'S8'><span>Solve the FBA problem with added constraints C*v &gt;= d , with or without objective function</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span style="color: rgb(2, 128, 9);">%modelConstrainedAb = changeObjective(modelConstrainedAb, 'DM_atp_c_');</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >FBAsolution = optimizeCbModel(modelConstrainedAb,</span><span style="color: rgb(170, 4, 249);">'max'</span><span >,1e-6);</span></span></div></div></div><div  class = 'S8'><span>Check the values of the added fluxes</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span >FBAsolution.x(rxnInd)</span></span></div><div  class = 'S7'><div class="inlineElement eoOutputWrapper embeddedOutputsMatrixElement" uid="80522B2E" data-testid="output_8" data-width="428" style="width: 458px;"><div class="matrixElement veSpecifier"><div class="veVariableName variableNameElement" style="width: 428px;"><span>ans = </span><span class="veVariableValueSummary"></span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:72,&quot;totalColumns&quot;:&quot;1&quot;,&quot;totalRows&quot;:&quot;3&quot;,&quot;charsPerColumn&quot;:10}"><div class="variableValue" style="width: 74px;">    2.0508
    0.2334
    0.3898
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div></div></div></div><div  class = 'S8'><span>Therefore, when you solve the FBA problem with this last constraint, the sum of flux values associated with these three reactions should be greater than the value of </span><span style=' font-style: italic;'>d</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S6'><span style="white-space: normal"><span >sum(c*FBAsolution.x(rxnInd))</span></span></div><div  class = 'S7'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>ans = 2.6740</div></div></div><div class="inlineWrapper"><div  class = 'S16'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">return</span><span >;</span></span></div></div></div><h2  class = 'S4'><span>CRITICAL STEP: Collection of data and conversion of experimental fluxes (Timing: 4-8 weeks)</span></h2><div  class = 'S3'><span>The most time consuming step when imposing constraints is the collection of required information. Depending on the available experimental literature, it can take between 4-8 weeks to retrieve the biomass composition and the turnover rates of the different biomass precursors. It is crucial to correctly convert the obtained data into the corresponding fluxes. It is recommended to first define the flux unit you wish to use. A common unit used for prokaryotic models is micromol per gramDryWeight per hour (</span><span texencoding="\mu mol/gDW/hr" style="vertical-align:-5px"><img src="" width="95.5" height="19" /></span><span>). However, in the experimental literature, a wide range of units is provided. Therefore, after each conversion, it is strongly recommended to double check the calculations to avoid modelling artifacts. Once all the constraints are available, it can take less than 5 minutes to impose the constraints on the corresponding reaction bounds, according to the information provided in this tutorial. </span></div><h2  class = 'S4'><span>ANTICIPATED RESULTS</span></h2><div  class = 'S3'><span>After imposing the above constraints, we can now test the likely outcome of an optimisation problem using a constraint-based model. For example, we can take advantage of sparseFBA to identify the minimal set of essential reactions required to fulfill a certain objective function (e.g. </span><a href = "http://www.vmh.life/#reaction/DM_atp_c_"><span>DM_atp_c_</span></a><span>).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >originalTest = model; </span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >originalTest = changeObjective(originalTest , </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >[vSparseOriginal, sparseRxnBoolOriginal, essentialRxnBoolOriginal]  = sparseFBA(originalTest);</span></span></div></div></div><div  class = 'S8'><span>Dis</span><span>play the number of essential reactions thata re required to carry flux to fulfill the objective function:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >cnt=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >n=1:length(originalTest.rxns)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >essentialRxnBoolOriginal(n,1)==true</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >        cnt=cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span > fprintf(</span><span style="color: rgb(170, 4, 249);">'%g%s\n'</span><span >,cnt,</span><span style="color: rgb(170, 4, 249);">' reactions essential to fulfill the objective function DM_atp_c_'</span><span >)</span></span></div></div></div><div  class = 'S8'><span>In the absence of constraints, the minimal set of reactions required to maximise the objective function is 111 essential reactions.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >constrainedTest = modelConstrainedAb; </span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >constrainedTest = changeObjective(constrainedTest, </span><span style="color: rgb(170, 4, 249);">'DM_atp_c_'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span >[vSparseConstrained, sparseRxnBoolConstrained, essentialRxnBoolConstrained]  = sparseFBA(constrainedTest);</span></span></div></div></div><div  class = 'S8'><span>Dis</span><span>play the number of essential reactions thata re required to carry flux to fulfill the objective function:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S9'><span style="white-space: normal"><span >cnt=0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">for </span><span >n=1:length(constrainedTest.rxns)</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">if </span><span >essentialRxnBoolConstrained(n,1)==true</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >        cnt=cnt+1;</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span >    </span><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S10'><span style="white-space: normal"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: normal"><span > fprintf(</span><span style="color: rgb(170, 4, 249);">'%g%s\n'</span><span >,cnt,</span><span style="color: rgb(170, 4, 249);">' reactions essential to fulfill the objective function DM_atp_c_'</span><span >)</span></span></div></div></div><div  class = 'S8'><span>After the addition of constraints, the minimal set of reactions required is increased to 172 essential reactions. Therefore, in this example, it is useful to integrate cell-type specific constraints to further define the minimal set of essential reactions. In most cases, constraints also allow us to alter the feasible solution space to obtain fluxes that better agree with the known physiology of the cell type. </span></div><h2  class = 'S4'><span>REFERENCES</span></h2><div  class = 'S3'><span>1. Feist, A.M. and Palsson, B.Ø. The biomass objective function. Current Opinion in Microbiology. 13(3), 344–349 (2010).  </span></div><div  class = 'S3'><span>2. Kuhar, M.J. On the use of protein turnover and half-lives. Neuropsychopharmacology. 34(5), 1172–1173 (2008). </span></div><div  class = 'S3'><span>3. Sokoloff, L. et al. The [14C]deoxyglucose method for the measurement of local cerebral  glucose utilization: theory, procedure, and normal values in the  conscious and anesthetized albino rat. J Neurochem. 28(5):897-916 (1977).</span></div><div  class = 'S3'><span>4. Thiele, I.  and Palsson B.Ø. A protocol for generating a high-quality genome-scale metabolic reconstruction. Nat. Protocols. 5(1), 93–121(2010).</span></div><div  class = 'S3'><span>5. Zhang, J. and Liu, Q. Cholesterol metabolism and homeostasis in the brain. Protein Cell. 6(4), 254-64 (2015). </span></div><div  class = 'S3'><span>6. Martinez-Vicente, M. Neuronal mitophagy in neurodegenerative diseases. Front. Mol. Neurosci. 8, 10:64 (2017).</span></div><div  class = 'S3'><span>7. Schulze, H. et al. Principles of lysosomal membrane degradation: cellular topology and biochemistry of lysosomal lipid degradation. Biochim. Biophys. Acta. 1793(4), 674-83 (2009). </span></div><div  class = 'S3'><span>8. Lajtha, A. and Sylvester, V. Handbook of Neurochemistry and Molecular Neurobiology. Springer; 2008. Available </span><a href = "http://www.springer.com/de/book/9780387303512"><span>here</span></a><span>. </span></div><div  class = 'S3'></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% Adding biological constraints to a flux balance model
% *Note: This tutorial is a draft and needs completion. Contributions welcome!*
%% Authors: Diana C. El Assal and Ronan M.T. Fleming,  Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Luxembourg.
%% Reviewers: 
% Anne Richelle, Systems Biology and Cell engineering, University of California 
% San Diego
% 
% Almut Heinken, Molecular Systems Physiology Group, LCSB, University of Luxembourg
%% INTRODUCTION
% A metabolic model can be converted into a condition-specific model based on 
% the imposition of experimentally derived constraints. Constraints can be defined, 
% for example, by imposing upper and lower flux bounds for each reaction. There 
% are several types of constraints that can be imposed in a metabolic model and 
% that represent specific intra- and extracellular conditions, such as biomass 
% maintenance requirements, environmental constraints, or maximum enzyme capacities. 
% 
% In general, biomass constraints [1] are added as part of a biomass reaction. 
% In some instances, however, a cell-type (e.g. neurons) does not divide, but 
% is only required to turn over its biomass components. This tutorial is particularly 
% relevant for such cases. Turnover rates are commonly expressed as half-lives 
% (${t_{1/2}$) and represent the time required for half of the biomass precursor 
% to be replaced [2]. 
% 
% Using the experimental literature, metabolite ${t_{1/2}$ were collected and 
% converted into turnover rates ($\lambda$):
% 
% $$\text{(1)}\ \lambda=\frac{ln(2)}{t_{1/2}}\\$$
% 
% We consider a biochemical network of  m  molecular species and  n  biochemical 
% reactions. The biochemical network is mathematically represented by a stoichiometric 
% matrix $S\in\mathcal{Z}^{m\times n}$.  After calculating $\lambda$, we integrate 
% it into the steady-state equation:
% 
% $$\text{(2)}\ Sv=\frac{dx}{dt}\\$$
% 
% The steady-state flux vector $v$ and the change in abundance over time ($\frac{dx}{dt}\\$) 
% share the same units, with  $x$  being the abundance of the corresponding biomass 
% precursor. Equation (2) can thus be re-written: 
% 
% $$\text{(3)}\ Sv=\frac{dx}{dt}=\lambda x\\$$
%% PROCEDURE
% Initialize the Cobra Toolbox using the |initCobraToolbox| function.

initCobraToolbox(false) % false, as we don't want to update
%% *Setting the* optimization *solver*

changeCobraSolver('gurobi','LP');
%% 
% Here, we use Recon2.0 model (distributed by the toolbox) for illustration, 
% although any model can be used. 

modelFileName = 'Recon2.0model.mat';
modelDirectory = getDistributedModelFolder(modelFileName); %Look up the folder for the distributed Models.
modelFileName= [modelDirectory filesep modelFileName]; % Get the full path. Necessary to be sure, that the right model is loaded
model = readCbModel(modelFileName);
modelOrig = model;
%% *1. Environmental constraints*
% 
% 
% Environmental constraints are typically related to nutrient availability (e.g., 
% glucose and oxygen). They can be defined using the function _changeRxnBounds_ 
% to set the minimal and maximal uptake and/or secretion rates possible in a specific 
% condition. For example, in the caudate-putamen of the conscious rat, glucose 
% consumption rate was found to range between -12.00 and -11.58  $\mu mol/gDW/hr$[3]. 
% Therefore, the lower bound of the glucose exchange reaction (<http://www.vmh.life/#reaction/EX_glc(e) 
% EX_glc(e)>) can be set as follows:

modelConstrained = model;
modelConstrained.c = 0*modelConstrained.c; % remove any objective function
modelConstrained = changeRxnBounds(modelConstrained, 'EX_glc(e)', -12, 'l');
%% 
% Optionally, to further constrain the model, an upper bound can also be imposed 
% to force the model to take up between 11.58 and 12 units of glucose 

modelConstrained = changeRxnBounds(modelConstrained, 'EX_glc(e)', -11.58, 'u');
%% *2. Internal enzymatic constraints* 
% 
% 
% By convention, the bounds set on reaction rates in a metabolic model range 
% from -1000 to 1000 and from 0 to 1000 for reversible and irreversible reactions, 
% respectively [4]. Actually, the rate of a reaction is related to the activity 
% of the enzyme catalyzing this reaction.  Therefore, internal enzymatic constraints 
% can be used to define the maximum capacity of a specified enzyme to catalyze 
% a reaction (_Vmax_). For example, assuming that the reaction catalyzed by fructose-bisphosphate 
% aldolase (<http://www.vmh.life/#reaction/FBA FBA>) has a _Vmax_ of 128 units 
% in our specific cell type, we can then add this constraint on the corresponding 
% internal reaction FBA, as an upper bound. 

modelConstrained = changeRxnBounds(modelConstrained, 'FBA', 128, 'u');
%% 
% Optionally, if the reaction is reversible, the same constraint can be set 
% as the lower bound, but with opposite signs.
% 
% 

modelConstrained = changeRxnBounds(modelConstrained, 'FBA', -128, 'l');
%% 
% 
%% *3. Constraints associated with biomass*
% In general, biomass constraints [1] are added as part of a biomass reaction 
% by defining stoichiometric coefficients for each biomass precursor. For dividing 
% cell types, the generic human biomass reaction available in Recon2 is formulated 
% as follows: 

printRxnFormula(modelConstrained, 'biomass_reaction');
%% *3.1 Biomass reaction*
% 20.6508 h2o[c] + 20.7045 atp[c] + 0.38587 glu_L[c] + 0.35261 asp_L[c] + 0.036117 
% gtp[c] + 0.50563 ala_L[c] + 0.27942 asn_L[c] + 0.046571 cys_L[c] + 0.326 gln_L[c] 
% + 0.53889 gly[c] + 0.39253 ser_L[c] + 0.31269 thr_L[c] + 0.59211 lys_L[c] + 
% 0.35926 arg_L[c] + 0.15302 met_L[c] + 0.023315 pail_hs[c] + 0.039036 ctp[c] 
% + 0.15446 pchol_hs[c] + 0.055374 pe_hs[c] + 0.020401 chsterol[c] + 0.002914 
% pglyc_hs[c] + 0.011658 clpn_hs[c] + 0.009898 dgtp[n] + 0.009442 dctp[n] + 0.013183 
% datp[n] + 0.053446 utp[c] + 0.013091 dttp[n] + 0.27519 g6p[c] + 0.12641 his_L[c] 
% + 0.15967 tyr_L[c] + 0.28608 ile_L[c] + 0.54554 leu_L[c] + 0.013306 trp_L[c] 
% + 0.25947 phe_L[c] + 0.41248 pro_L[c] + 0.005829 ps_hs[c] + 0.017486 sphmyln_hs[c] 
% + 0.35261 val_L[c]  -> 20.6508 h[c] + 20.6508 adp[c] + 20.6508 pi[c] 
% 
% Any changes or adaptations can be introduced by adding a new formulation of 
% the biomass function, using the function _addReaction_. For example, one can 
% add the following new biomass reaction named _biomassReactionLipids_:

modelConstrained = addReaction(modelConstrained, 'biomasReactionLipids',  '20.6508 h2o[c] + 20.7045 atp[c] + 0.15446 pchol_hs[c] + 0.055374 pe_hs[c] + 0.020401 chsterol[c] + 0.011658 clpn_hs[c] + 0.005829 ps_hs[c] + 0.017486 sphmyln_hs[c] -> 20.6508 h[c] + 20.6508 adp[c] + 20.6508 pi[c]');
%% *4. Biomass maintenance constraints*
% To represent biomass maintenance (e.g. in neurons), the minimal biomass maintenance 
% requirements can be used to set the corresponding constraints. Using the neurobiochemical 
% literature, the degradation pathways for each biomass precursor has to be identified 
% and the corresponding first reactions of these degradation pathways need to 
% be mapped to ReconX. Using the fractional composition and the turnover rate 
% of each biomass precursor, corresponding reaction rates ($\mu mol/gDW/hr$) are 
% calculated as described above and in Table 1. These reaction rates represent 
% the minimal requirements for biomass maintenance of neurons in the human grey 
% matter and must therefore be imposed as a lower bound on the corresponding degradation 
% reaction(s) of the different lipids, amino acids, and nucleic acids. 
% 
% <#null >
% 
% *Table 1: The minimum metabolic maintenance requirement for neurons.* This 
% is a coarse-grained approximation of neuronal lipid, amino acid, and nucleic 
% acid maintenance requirements converted into $\mu mol/gDW/hr$.  
% 
% 
% 
% **cardiolipin is also known as diphosphatidylglycerol
% 
% _*Calculation example of the minimal cholesterol maintenance requirement*_ 
% 
% _*i. Identify metabolite abundance*_
% 
% Assuming that the specific tissue type has a total dry weight lipid composition 
% of 39.6%. This means that there is 0.396gLipid/gDWtissue. If cholesterol has 
% a molar composition of 31.3%, then in total there is 0.124g cholesterol per 
% gDW of tissue.
% 
% $$\frac{31.3*39.6*1gDW}{100*100}=0.124gDW\\\$$

abundance = (31.3*39.6*1)/(100*100);
%% 
% _*ii. Calculate the molar abundance*_
% 
% In the experimental literature, cholesterol was also found to have a molar 
% mass ($M$) of 386$g/mol$. Using equation (4), we can now convert the abundance 
% ($m$) into molar units ($n$). 
% 
% $$\text{(4)}\ n=\frac{m}{M}$$

M = 386; %g/mol
n = (abundance*1000000)/M; %micromol
%% 
% _*iii. Calculate the corresponding flux value*_
% 
% Finally, we know that in the brain, cholesterol has a very slow turnover and 
% a ${t_{1/2}$ of 4320 hours. Using equation (3), we can now calculate the minimal 
% cholesterol maintenance requirement in flux units ($v_1$). 

halfLife = 4320;
turnover = log(2)/halfLife;
v1 = n * turnover
%% 
% The minimal cholesterol maintenance requirement was calculated to be 0.0515$\mu 
% mol/gDW/hr$(Table 1). This value can now be used as a lower bound in the corresponding 
% reaction.
%% *4.1.* Identification of degradation reactions for a biomass maintenance precursor
% As previously mentioned, the degradation pathways for each biomass precursor 
% can be identified using literature and the minimal maintenance requirements 
% defined in section 4. can be used to constrain the first reactions of these 
% degradation pathways. However, the identification of these reactions and the 
% set-up of the associated constraints is not always straightforward. The following 
% section presents the different common cases that can be encountered.
%% A. Single irreversible degradation reaction
% In cases where only a single irreversible degradation reaction exists for 
% a biomass maintenance precursor, the imposition of the constraint is straightforward. 
% For example, the major cholesterol excretion pathway in the brain involves the 
% hydroxylation of cholesterol into the oxysterol 24-hydroxycholesterol. Only 
% a subset of neurons express this 24-hydroxylase enzyme (<http://www.vmh.life/#human/all/P45046A1r 
% P45046A1r>) and it is mainly found in dendrites and somata, rather than in axons 
% or presynaptic terminals (reviewed in [5]). 
% 
% Therefore, add a lower bound ($v_1$, calculated above) on P45046A1r (Table 
% 1). 

modelConstrained = changeRxnBounds(modelConstrained, 'P45046A1r', v1, 'l');
%% B. Single degradation reaction does not exist biochemically
% A degradation reaction might not exist for a given biomass maintenance precursor. 
% For example, the phospholipid cardiolipin is mainly present in the inner mitochondrial 
% membrane, where it regulates the stability of the mitochondrial membrane protein 
% complexes [6]. As part of mitochondria, cardiolipin reaches the lysosome during 
% macroautophagy (reviewed in [7]). It is then degraded to form the negatively 
% charged bis(monoacylglycero)phosphate (BMP) on internal membranes.  
% 
% If the corresponding demand reaction does not exist in the model, a demand 
% reaction can be added using the function _addDemandReaction_. 

modelConstrained = addDemandReaction(modelConstrained, 'clpn_hs[c]'); 
%% 
% Now, the constraint for cardiolipin (Table 1) can be imposed on the corresponding 
% demand reaction

modelConstrained = changeRxnBounds(modelConstrained, 'DM_clpn_hs[c]', 0.001, 'l');
%% C. Single reversible degradation reaction
% In cases where the biomass precursor degradation reaction is reversible, it 
% first needs to be split into two irreversible reactions. To this end, define 
% the set of reversible degradation reactions (sRxns) and convert them into two 
% irreversible reactions (i.e., _sRXN_b_ and _sRXN_f_, respectively backward and 
% forward reactions) using _convertToIrreversible:_
% 
% *i. Split the reversible reactions into irreversible*
% 
% Select a set of reversible degradation reactions 

sRxns = {'ASPTA'; 'GHMT2r'};
%% 
% Copy the original model, except split a specific list of reversible degradation 
% reactions into irreversible. 
% 
% Split sRxns into irreversible reactions

[modelIrrev] = convertToIrreversible(modelConstrained,'sRxns',sRxns);
%% 
% You can check if the conversion as been done properly by searching the split 
% reactions

for j=1:length(sRxns)
    if isempty(findRxnIDs(modelIrrev, [sRxns{j} '_f']))
        error('Forward reaction not found')
    end
    if isempty(findRxnIDs(modelIrrev, [sRxns{j} '_b']))
        error('Reverse reaction not found')
    end
end
%% 
% *ii. Impose the calculated constraints*
% 
% Examples are given for aspartate and glutamate (Table 1). 

constraints = [1.590; 1.146];
%% 
% You can also identify the new reaction names as follows:

rxns = setdiff(modelIrrev.rxns, modelConstrained.rxns)
%% 
% Using this list (rxns), manually identify the corresponding reactions that 
% should be constrained

splitRxns= {'ASPTA_f'; 'GHMT2r_f'};
%% 
% Identify the indices of these split reactions in the new model, using _findRxnIDs_

ind = findRxnIDs(modelIrrev, splitRxns);
%% 
% Now, impose the constraints using a for loop. Note that you can easily account 
% for experimental errors by defining a percentage error (e.g., _expError = 0.25_) 
% for the constraint values. 

expError = 0.25;
modelConstrained = modelIrrev;
for i = 1:length(splitRxns)
    modelConstrained = changeRxnBounds(modelConstrained, splitRxns{i,1},...
        constraints(i,1)-constraints(i,1)*expError, 'l');
end
%% D. Multiple irreversible degradation reactions
% In some cases, several degradation pathways may be available for one biomass 
% precursor. For example, in the brain, phosphatidylcholine (PC) can be degraded 
% by 3 different metabolic pathways [8]:
%% 
% * <http://www.vmh.life/#reaction/PCHOLP_hs PCHOLP_hs>: Phospholipase D acts 
% on the choline/phosphate bond of PC to form choline and phosphatidic acid.
% * <http://www.vmh.life/#reaction/PLA2_2 PLA2_2>: Phospholipase A2 acts on 
% the bond between the fatty acid and the hydroxyl group of PC to form a fatty 
% acid (e.g. arachidonic acid or docosahexaenoic acid) and lysophosphatidylcholine. 
% * <http://www.vmh.life/#reaction/SMS SMS>: Ceramide and PC can also be converted 
% to sphingomyelin by sphingomyelin synthetase. 
%% 
% Define the set of potential reactions associated with the degradation of PC

multipleRxnList={'PCHOLP_hs', 'PLA2_2', 'SMS'};
%% 
% Make sure that all the reactions are irreversible. The lower bounds should 
% be 0 and the upper bounds 1000.  

modelConstrained.lb(findRxnIDs(modelConstrained, multipleRxnList));
%% 
% ans =
% 
% 0
% 
% 0
% 
% 0

modelConstrained.ub(findRxnIDs(modelConstrained, multipleRxnList));
%% 
% ans =
% 
% 1000
% 
% 1000
% 
% 1000
% 
% Constrain the weighted sum of fluxes to be above a lower bound (e.g. value 
% of the maintenance requirement of PC in Table 1: d = 2.674 $\mu mol/gDW/hr$). 
% The weight for each reaction are defined in _c_.

c=[1,1,1];
d=2.674; 
ineqSense='G';
modelConstrainedAb=constrainRxnListAboveBound(modelConstrained,multipleRxnList,c,d,ineqSense);
rxnInd = findRxnIDs(modelConstrainedAb, multipleRxnList);
%% 
% Check the constraints are there

[nMet,nRxn]=size(modelConstrainedAb.S);
modelConstrainedAb
%% 
% Solve the FBA problem with added constraints C*v >= d , with or without objective 
% function

%modelConstrainedAb = changeObjective(modelConstrainedAb, 'DM_atp_c_');
FBAsolution = optimizeCbModel(modelConstrainedAb,'max',1e-6);
%% 
% Check the values of the added fluxes

FBAsolution.x(rxnInd)
%% 
% Therefore, when you solve the FBA problem with this last constraint, the sum 
% of flux values associated with these three reactions should be greater than 
% the value of _d_

sum(c*FBAsolution.x(rxnInd))
return;
%% CRITICAL STEP: Collection of data and conversion of experimental fluxes (Timing: 4-8 weeks)
% The most time consuming step when imposing constraints is the collection of 
% required information. Depending on the available experimental literature, it 
% can take between 4-8 weeks to retrieve the biomass composition and the turnover 
% rates of the different biomass precursors. It is crucial to correctly convert 
% the obtained data into the corresponding fluxes. It is recommended to first 
% define the flux unit you wish to use. A common unit used for prokaryotic models 
% is micromol per gramDryWeight per hour ($\mu mol/gDW/hr$). However, in the experimental 
% literature, a wide range of units is provided. Therefore, after each conversion, 
% it is strongly recommended to double check the calculations to avoid modelling 
% artifacts. Once all the constraints are available, it can take less than 5 minutes 
% to impose the constraints on the corresponding reaction bounds, according to 
% the information provided in this tutorial. 
%% ANTICIPATED RESULTS
% After imposing the above constraints, we can now test the likely outcome of 
% an optimisation problem using a constraint-based model. For example, we can 
% take advantage of sparseFBA to identify the minimal set of essential reactions 
% required to fulfill a certain objective function (e.g. <http://www.vmh.life/#reaction/DM_atp_c_ 
% DM_atp_c_>).

originalTest = model; 
originalTest = changeObjective(originalTest , 'DM_atp_c_');
[vSparseOriginal, sparseRxnBoolOriginal, essentialRxnBoolOriginal]  = sparseFBA(originalTest);
%% 
% Display the number of essential reactions thata re required to carry flux 
% to fulfill the objective function:

cnt=0;
for n=1:length(originalTest.rxns)
    if essentialRxnBoolOriginal(n,1)==true
        cnt=cnt+1;
    end
end
 fprintf('%g%s\n',cnt,' reactions essential to fulfill the objective function DM_atp_c_')
%% 
% In the absence of constraints, the minimal set of reactions required to maximise 
% the objective function is 111 essential reactions.

constrainedTest = modelConstrainedAb; 
constrainedTest = changeObjective(constrainedTest, 'DM_atp_c_');
[vSparseConstrained, sparseRxnBoolConstrained, essentialRxnBoolConstrained]  = sparseFBA(constrainedTest);
%% 
% Display the number of essential reactions thata re required to carry flux 
% to fulfill the objective function:

cnt=0;
for n=1:length(constrainedTest.rxns)
    if essentialRxnBoolConstrained(n,1)==true
        cnt=cnt+1;
    end
end
 fprintf('%g%s\n',cnt,' reactions essential to fulfill the objective function DM_atp_c_')
%% 
% After the addition of constraints, the minimal set of reactions required is 
% increased to 172 essential reactions. Therefore, in this example, it is useful 
% to integrate cell-type specific constraints to further define the minimal set 
% of essential reactions. In most cases, constraints also allow us to alter the 
% feasible solution space to obtain fluxes that better agree with the known physiology 
% of the cell type. 
%% REFERENCES
% 1. Feist, A.M. and Palsson, B.Ø. The biomass objective function. Current Opinion 
% in Microbiology. 13(3), 344–349 (2010).  
% 
% 2. Kuhar, M.J. On the use of protein turnover and half-lives. Neuropsychopharmacology. 
% 34(5), 1172–1173 (2008). 
% 
% 3. Sokoloff, L. et al. The [14C]deoxyglucose method for the measurement of 
% local cerebral  glucose utilization: theory, procedure, and normal values in 
% the  conscious and anesthetized albino rat. J Neurochem. 28(5):897-916 (1977).
% 
% 4. Thiele, I.  and Palsson B.Ø. A protocol for generating a high-quality genome-scale 
% metabolic reconstruction. Nat. Protocols. 5(1), 93–121(2010).
% 
% 5. Zhang, J. and Liu, Q. Cholesterol metabolism and homeostasis in the brain. 
% Protein Cell. 6(4), 254-64 (2015). 
% 
% 6. Martinez-Vicente, M. Neuronal mitophagy in neurodegenerative diseases. 
% Front. Mol. Neurosci. 8, 10:64 (2017).
% 
% 7. Schulze, H. et al. Principles of lysosomal membrane degradation: cellular 
% topology and biochemistry of lysosomal lipid degradation. Biochim. Biophys. 
% Acta. 1793(4), 674-83 (2009). 
% 
% 8. Lajtha, A. and Sylvester, V. Handbook of Neurochemistry and Molecular Neurobiology. 
% Springer; 2008. Available <http://www.springer.com/de/book/9780387303512 here>. 
% 
%
##### SOURCE END #####
-->
</div></body></html>